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The NP-hard problem of optimizing a quadratic form over the unimodular vector set arises in radar 
code design scenarios as well as other active sensing and communication applications. To tackle this 
problem (which we call unimodular quadratic programming (UQP)), several computational approaches 
■ are devised and studied. A specialized local optimization scheme for UQP is introduced and shown to 

yield superior results compared to general local optimization methods. Furthermore, a monotonically 



CZ2 

error-bound improving technique (MERIT) is proposed to obtain the global optimum or a local optimum 
of UQP with good sub-optimality guarantees. The provided sub-optimality guarantees are case-dependent 
and generally outperform the n/4 approximation guarantee of semi-definite relaxation. Several numerical 
I/") ' examples are presented to illustrate the performance of the proposed method. The examples show that 

for cases including several matrix structures used in radar code design, MERIT can solve UQP efficiently 
CO . in the sense of sub-optimality guarantee and computational time. 

o 

Index Terms 

' radar codes, unimodular codes, quadratic programming. 

I. Introduction 

Unimodular codes are used in many active sensing and communication systems mainly as a result of 
the their optimal (i.e. unity) peak-to-average -power ratio (PAR). The design of such codes can be often 
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formulated as the optimization of a quadratic form (see sub-section II-AI for examples). Therefore, we 
will study the problem 

UQP: max s H Rs (1) 

sefi n 

where R £ C nxn is a given Hermitian matrix, (.) H denotes the vector/matrix Hermitian transpose, £1 
represents the unit circle, i.e. Q = {s € C : \s\ = 1} and UQP stands for Unimodular Quadratic 
Program(ming). 

A. Motivating Applications 

To motivate the UQP formulation considered above, we present four scenarios in which a design 
problem in active sensing or communication boils down to an UQP. 

• Designing codes that optimize the SNR or the CRLB: We consider a monostatic radar which transmits 
a linearly encoded burst of pulses. The observed backscattered signal v can be written as (see, e.g. CQ): 

v = a(cQp)+w, (2) 

where a represents both channel propagation and backscattering effects, w is the disturbance/noise com- 
ponent, c is the unimodular vector containing the code elements, p = (1, e i 27r /rf T »- ) . . . ^ e j2n(n-i)f d T r ^T 
is the temporal steering vector with and T r being the target Doppler frequency and pulse repetition 
time, respectively, and the symbol stands for the Hadamard (element- wise) product of matrices. 

Under the assumption that id is a zero-mean complex-valued circular Gaussian vector with known 
positive definite co variance matrix K[ww H ] = M, the signal-to-noise ratio (SNR) is given by [|2l 

SNR = \a\ 2 c H Rc (3) 

where R = M" 1 (pp H )* with (.)* denoting the vector/matrix complex conjugate. Therefore, the 
problem of designing codes optimizing the SNR of the radar system can be formulated directly as an 
UQP. Additionally, the Cramer-Rao lower bound (CRLB) for the target Doppler frequency estimation 
(which yields a lower bound on the variance of any unbiased target Doppler frequency estimator) is given 
by m 

CRLB = (2\a\ 2 (cQpQu) H M' 1 (cQpQu)y 1 (4) 
= (2\a\ 2 c H tf C y l 

where u = (0, j2u:T r , ■ ■ ■ ,j2vr(n - l)T r ) T and R' = M' 1 {pp H f (uu H )*. Therefore the 
minimization of CRLB can also be formulated as an UQP. For the simultaneous optimization of SNR 
and CRLB see (2l. 
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• Synthesizing cross ambiguity functions (CAFs): The ambiguity function (which is widely used in 
active sensing applications 011) represents the two-dimensional response of the matched filter to a 
signal with time delay r and Doppler frequency shift /. The more general concept of cross ambiguity 
function occurs when the match filter is replaced by a mismatched filter. The cross ambiguity function 
(CAF) is defined as 

POO 

X (tJ)= / u(t)v*{t + T)e^f l dt (5) 



where u(t) and v(t) are the transmit signal and the receiver filter, respectively (the ambiguity function 
is obtained from © with v(t) = u(t)). In several applications u(t) and v(t) are given by: 



u(t) = v(t) = ^2ykPk(t) (6) 

k=l k=l 

where {pk{t)} are pulse-shaping functions (with the rectangular pulse as a common example), and 

x = (xi ■ ■ -Xnf, y = {yi---y n ) T (7) 

are the code and, respectively, the filter vectors. The design problem of synthesizing a desired CAF has a 
small number of free variables (i.e. the entries of the vectors x and y) compared to the large number of 
constraints arising from two-dimensional matching criteria (to a given |x( r > /)!)• Therefore, the problem 
is generally considered to be difficult and there are not many methods to synthesize a desired (cross) 
ambiguity function. Below, we describe briefly the cyclic approach of @ for CAF design. 

The problem of matching a desired |x( t j /)| = d(r, f) can be formulated as the minimization of the 
criterion 

poo roo 2 



drdf (8) 



/OO POO 
/ w(r,f) d(r,f)eP^ -y H J(T,f)x 
-oo J —oo 

where «7(r, /) G C nxn is given, w(t, f) is a weighting function that specifies the CAF area which needs 
to be emphasized and <f>(r, f) represent auxiliary phase variables. It is not difficult to see that for fixed 
x and y, the minimizer <p(r,f) is given by (j)(r,f) = &ig{y H J(r, f)x}. For fixed 0(r, /) and x, the 
criterion g can be written as 

g(y) = y H D lV - y H B H x - x 11 By + consti (9) 

= (y - DY 1 B H x) H D 1 (y - D^ 1 B H x) +const 2 

where B and D\ are given matrices in C nxn 0. Due to practical considerations, the transmit coefficients 
{xk} must have low PAR values. However, the receiver coefficients {yk} need not be constrained in such 
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a way. Therefore, the minimizer y of g(y) is given by y = D x l B H x. Similarly, for fixed 4>(t, /) and 
y, the criterion g can be written as 

g(x) = x H D 2 x - x H By - y H B H x + const 3 (10) 

where D 2 € C nxn is given 0. If a unimodular code vector x is desired then the optimization of g(x) 
is an UQP as g{x) can be written as 

D 2 -By \ I e>f>x 
g(x) = . +const 3 (11) 




-{By) H n / \ p3V 
where ip G [0, 2ir) is a free phase variable. 

• Steering vector estimation in adaptive beamforming: Consider a linear array with n antennas. The 
output of the array at time instant k can be expressed as ||6l 

x k = s k a + n k (12) 

with {s k } being the signal waveform, a the associated steering vector (with \[a]i\ = 1, 1 < I < ri), and 
n k the vector accounting for all independent interferences. 

The true steering vector is usually unknown in practice, and it can therefore be considered as an uni- 
modular vector to be determined |]7]. Define the sample covariance matrix of {x k } as R = ^ Y2k=i x kX k 
where T is the number of training data samples. Assuming some prior knowledge on a (which can be 
represented by arg(a) being in a given sector 0), the problem of estimating the steering vector can be 
formulated as [8| 



min a H R a (13) 

a 

s.t. arg(a) G G, 

hence an UQP-type problem. Such problems can be tackled using general local optimization techniques 
or the optimization scheme introduced in Section In] 

• Maximum likelihood (ML) detection of unimodular codes: Assume the linear model 

y = Qs + n (14) 

where Q represents a multiple-input multiple-output (MIMO) channel, y is the received signal, n is the 
additive white Gaussian noise and s contains the unimodular symbols which are to be estimated. The 
ML detection of s may be stated as 

s M l = arg min \\y - Qs\\ 2 (15) 

sen™ 
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It is straightforward to verify that the above optimization problem is equivalent to the UQP Q: 



where 



min s H Rs (16) 



Q"Q -Q"y 





( e>fa \ 


)•-! 






i & J 



-y H Q o 

and ip G [0, 27r) is a free phase variable. 



B. Related Work 

In |[T0l , the NP-hardness of UQP is proven by employing a reduction from an NP-complete matrix 
partitioning problem. The UQP in CO is often studied along with the following (still NP-hard) related 
problem in which the decision variables are discrete: 

m-UQP: m&xs H Rs (18) 

where Q m = {l,e?%,--- ,e^ ( - m ^ 1 ^}. Note that the latter problem coincides with the UQP in (Q} 
as m — > oo. The authors of ifTTI show that when the matrix R is rank-deficient (more precisely, when 
d =rank(i?) behaves like 0(1) with respect to the problem dimension) the m-UQP problem can be solved 
in polynomial-time and they propose a 0((?nn/2) 2d )-complexity algorithm to solve ([T8T l. However, such 
algorithms are not applicable to the UQP which corresponds to an infinite m. 

Studies on polynomial-time algorithms for UQP (and m-UQP) have been extensive (e.g. see [9]-[19] 
and the references therein). In particular, the semi-definite relaxation (SDR) technique has been one of 
the most appealing approaches to the researchers. To derive an SDR, we note that s H Rs = ti(s H Rs) = 
tr(Rss H ). Hence, the UQP can be rewritten as 

max tr(RS) (19) 
s 

s.t. s = ss 11 , s g n n . 

If we relax (fT9l by removing the rank constraint on S and the unimodularity constraint on s then the 
result is a semi-definite program: 

SDP: maxtr( J R5) (20) 
s.t. [S]k,k = 1, 1 < k < n, 
S is positive semi-definite. 
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The above SDP can be solved in polynomial time using interior-point methods |[T31 . The approximation 
of the UQP solution based on the SDP solution can be accomplished in several ways. For example, 
we can approximate the phase values of the solution s using a rank-one approximation of S. A more 
effective approach for guessing s is based on randomized approximations (see iTTOl . |[T6l and ifTTl ). A 
detailed guideline for randomized approximation of the UQP solution can be found in fTTH . In addition, 
we refer the interested reader to the survey of the rich literature on SDR in lfl~8l . 

Analytical assessments of the quality of the UQP solutions obtained by SDR and randomized approx- 
imation are available. Let vsdr be the expected value of the UQP objective at the obtained randomized 
solution. Let v op t represent the optimal value of the UQP objective. We have 

JV op t < V S DR < V opt (21) 

with the sub-optimality guarantee coefficient 7 = 7r/4 |[T0lllfT9ll . Note that the sub-optimality coefficient 
of the solution obtained by SDR can be arbitrarily close to ir/4 (e.g., see lfl~9l ). 

C. Contributions of this Work 

Besides SDR, the literature does not offer many other numerical approaches to tackle UQP. In this paper, 
a specialized local optimization scheme for UQP is proposed. The proposed computationally efficient 
local optimization approach can be used to tackle UQP as well as improve upon the solutions obtained 
by other methods such as SDR. Furthermore, a monotonically error-bound improving technique (called 
MERIT) is introduced to obtain the global optimum or a local optimum of UQP with good sub-optimality 
guarantees. Note that: 

• MERIT provides case-dependent sub-optimality guarantees. To the best of our knowledge, such 
guarantees for UQP were not known prior to this work. Using the proposed method one can generally 
obtain better performance guarantees compared to the analytical worst-case guarantees (such as 
7 = tt/4 for SDR). 

• The provided case-dependent sub-optimality guarantees are of practical importance in decision mak- 
ing scenarios. For instance in some cases the UQP solution obtained by SDR (or other optimization 
methods) might achieve good objective values. However, unless the goodness of the obtained solution 
is known (this goodness can be determined using the proposed bounds), the solution cannot be trusted. 

• Using MERIT, numerical evidence is provided to show that several UQPs (particularly those which 
occur in active sensing code design) can be solved efficiently without sacrificing the solution 
accuracy. 
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Finally, we believe that the general ideas of this work can be adopted to tackle m-UQP as the finite 
alphabet case of UQP. However, a detailed study of m-UQP is beyond the scope of this paper. 

D. Organization of the Paper 

The rest of this work is organized as follows. Section |ITJ discusses several properties of UQP. Section 
Hill introduces a specialized local optimization method. Section [TV] presents a cone approximation that is 
used in Section [V] to derive the algorithmic form of MERIT for UQP. Several numerical examples are 
provided in section [VlJ Finally, Section IVIII concludes the paper. 

Notation: We use bold lowercase letters for vectors/sequences and bold uppercase letters for matrices. 
(.) T denotes the vector/matrix transpose. 1 and are the all-one and all-zero vectors/matrices, e k is 
the k th standard basis vector in C n . \\x\\ n or the / n -norm of the vector x is defined as (£\ \x{k)\ n )^ 

where {x(k)} are the entries of x. The Frobenius norm of a matrix X (denoted by with entries 

i 

{X(k,l)} is equal to (^2 k ; \X (k, l)\ 2 ^j 2 . We use 9ft (X) to denote the matrix obtained by collecting 
the real parts of the entries of X. The matrix e? x is defined element-wisely as [e jX ] fci = e^ x ^ k - 1 . 
arg(.) denotes the phase angle (in radians) of the vector/matrix argument. E[.] stands for the expectation 
operator. Diag(.) denotes the diagonal matrix formed by the entries of the vector argument, whereas 
diag(.) denotes the vector formed by collecting the diagonal entties of the matrix argument. e>fc(X) 
represents the k th maximal eigenvalue of X. Finally, M. and C represent the set of real and complex 
numbers, respectively. 

II. Some Properties of UQP 

In this section, we study several properties of UQP. The discussed properties lay the grounds for a 
better understanding of UQP as well as the tools proposed to tackle it in the following sections. 

A. Basic Properties 

The UQP formulation in £[]) covers both maximization and minimization of quadratic forms (one can 
obtain the minimization of the quadratic form in CD) by considering — R in lieu of R). In addition, 
without loss of generality, the Hermitian matrix R can be assumed to be positive (semi)definite. If R is 
not positive (semi)definite, we can make it so using the diagonal loading technique (i.e. R <— R + XI 
where A > —a n (R)). Note that such a diagonal loading does not change the solution of UQP as 
s H (R + \I)s = s H Rs + \n. Next, we note that if s is a solution to UQP then e^s (for any £ [0, 2ir)) 
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is also a valid solution. To establish connections among different UQPs, Theorem 1 presents a bijection 
among the set of matrices leading to the same solution. 

Theorem 1. Let fC(s) represent the set of matrices Rfor which a given s £ £l n is the global optimizer 
of UQP Then 

1) IC(s) is a convex cone. 

2) For any two vectors Si,s 2 £ £l n , the one-to-one mapping (where s = s 2 ) 



It is interesting to note that in light of the above result, the characterization of the cone /C(s) for 
any given s = s leads to a complete characterization of all /C(s), s £ 0™, and thus solving any UQP. 
However, the NP-hardness of UQP suggests that such a characterization cannot be expected. Further 
discussions regarding the characterization of IC(s) are deferred to Section HVl 

B. Analytical Solutions to UQP 

There exist cases for which the analytical global optima of UQP are easy to obtain. In this sub-section, 
we consider two such cases which will be used later to provide an approximate characterization of /C(s). 
A special example is the case in which e J ^(-R) ( see the notation definition in II-Db is a rank-one matrix. 
More precisely, let R = R\ (ss ) where R\ is a real-valued Hermitian matrix with non-negative 
entties and s £ U n (a simple special case of this example is when R is a rank-one matrix itself). In 
this case, it can be easily verified that Ri £ /C(l nx i). Therefore, using Theorem 1 one concludes that 
R £ lC(s) i.e. s = s yields the global optimum of UQP. As another example, Theorem 2 considers the 
case for which several largest eigenvalues of the matrix R are identical. 

Theorem 2. Let R be a Hermitian matrix with eigenvalue decomposition R = UHU H . Suppose X is 
of the form 




(22) 



holds among the matrices in fC(si) and K,(s2). 



Proof: See the Appendix. 




O'n-m+l] ) 



(23) 



m times 



a\ > cr 2 > • • • > cr n _ m+ i 
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and let U m be the matrix made from the first m columns of U. Now suppose s £ f2 n lies in the linear 
space spanned by the columns of U m , i.e. there exists a vector a G C m such that 

s = U m a. (24) 

Then s is a global optimizer of UQP. 

Proof: Refer to the Appendix. ■ 

We end this section by noting that the solution to an UQP is not necessarily unique. For any set of 
unimodular vectors {si, S2, • • • , k < n, we can use the Gram-Schmidt process to obtain a unitary 
matrix U the first k columns of which span the same linear space as si, S2, • • • , Sfc- In this case, Theorem 
2 suggests a method to construct a matrix R (by choosing a £ with k identical largest eigenvalues) for 
which all si, S2, • • • , are global optimizers of the corresponding UQP. 

III. Specialized Local Optimization of UQP 

Due to its NP-hard nature, UQP has in general a highly multi-modal optimization objective. Finding 
and studying the local optima of UQP is not only useful to tackle the problem itself (particularly for 
UQP-related problems such as (fl"3T)), but also to improve the UQP approximate solutions obtained by 
SDR or other optimization techniques. In this section, we introduce a computationally efficient procedure 
to obtain a local optimum of UQP. 

Note that, while the risk for this to happen in practice is nearly zero, local optimization methods 
can in theory converge to a saddle point. Consequently, in the sequel we let L represent the set of all 
local optima and saddle points of UQP. Moreover, we assume that R is positive definite. Consider the 
following relaxed version of UQP: 

(RUQP) max K(sf#s 2 ) (25) 

We note that for fixed s 2 the maximizer of RUQP is given by 

ai =e jarg(K S2 )_ (26) 

Similarly, for any fixed si the maximizer of RUQP is given by 

S2 = e jarg(K Sl ). (27 ) 

In the following, we show that such a cyclic maximization of d25l ) can be used to find local optima of 
UQP. It is not difficult to see that the criterion in (|25T ) increases and is upper bounded (by ^ fe ; \ R(k, 
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through the iterations in (I26ll-(l271i. thus the said iterations are convergent in the sense of associated 
objective value. Next consider the identity 

25t(sf Rs 2 ) = s^Rsi + s%Rs 2 (28) 
- (ai - s 2 ) H R( Sl - a 2 ). 



Define e s = ||ai — s 2 \\ 2 an ^ suppose that s 2 is fixed and its associated optimal s\ is obtained by (126I ). 
It follows from ([28]) that 

K(af Rs 2 ) < 1 (af fiai + af fls 2 ) - ^a n (R). (29) 

Now suppose s 2 is the optimal vector in Q n obtained by (|2VT > for the above ai. Observe that ^(s^ 1 Rs 2 ) > 
K(«f and that K(sf iia' 2 ) > »(af Rs 2 ) > ^(sfRs 2 ) which imply 

5R(af i*a' 2 ) > - (af Bai + «^ifcj 2 ) (30) 



> R(afJk 2 ) + -^ ft (Ji). 



It follows from (f30]> that 



e 8 < —^\M(s?Rs' 2 )-$l(s?Rs 2 )\. (31) 
cr n (-K) 

The right-hand side of (OH vanishes through the cyclic minimization in (|26T)- (|27T ) which implies that e s 
converges to zero at the same time. Note that the above arguments can be repeated for fixed Si. We 
conclude that the iterations in (l26l>-(f27T> are convergent and also that they cannot converge to (si,s 2 ) 
with Si ^ s 2 . Moreover, as s H Rs = K(s^ Rs 2 ) for any Si = s 2 = s then any local optimum (ai,a 2 ) 
of RUQP satisfying si = s 2 = s yields a local optimum s of UQP. Based on the above discussions, the 
cyclic optimization of RUQP can be used to find local optima of UQP. Particularly, starting from any 
vector a(°) G Q n , the power method-like iterations 

g (t+l) =e iaxg( jRa «) (32) 

converge to an element in L. As an aside remark, we show that the objective of UQP is also increasing 
through the iterations of (|32l . Using (|28T ) with s\ = s^ t+1 \ and s 2 = (s( t+1 ) 7^ gW) implies that 

s (t + l)H Rs (t + l) > _ s (t)H Rs (t) +m{s{ t + l)H Rs (t) ) 

> s^ H RsV. (33) 

Note that while d32~l ) can obtain the local optima of UQP, it might not converge to every of them. To 
observe this, let s\ be a local optimum of UQP and initialize d32l) with s^> = s\. Let s 2 be another local 
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optimum of UQP but with a larger value of UQP than that at s{. Now one can observe from (l28l l that if £2 
is sufficiently close to si then the above power method-like iterations can move away from S\, meaning 
that they can converge to another local optimum of UQP with a larger value of the UQP objective than 
that at s[. Therefore, (l32l i bypasses some local optima of UQP with relatively small UQP objective values 
(which can be considered as an advantage compared to a general local optimization method). Moreover, 
one can note that there exist initializations for which (l32l i leads to the global optimum of UQP (i.e. the 
global optimum is not excluded from the local optima to which (l32l) can converge). 

Next, we observe that any s G L obtained by the above local optimization can be characterized by the 
equation 

arg(s) = aig(Rs). (34) 

We refer to the subset of L satisfying (l34l) as the hyper points of UQP. Note that if s € $7™ is a hyper 
point of UQP, then (l34b follows from the convergence of (l32l ). On the other hand, if (l34l is satisfied, 
it implies the convergence of the iterations in (|32l l and as a result s being a hyper point of UQP. The 
characterization given in (l34l is used below to motivate the characterization approach of Theorem 3. 

IV. Results on the cone /C(s) 

While a complete characterization of IC(s) cannot be expected (due to the NP-hardness of UQP), 
approximate characterizations of IC(s) axe possible. The goal of this section is to provide an approximate 
characterization of the cone /C(s) which can be used to tackle the UQP problem. Our main result is as 
follows: 

Theorem 3. For any given s = (e^ 1 , • • • ,e^ n ) T E Q n , let {B^i} be a set of matrices defined as 

B Ki = (e k ef + e ie * ) (ss H ) (35) 

and V s = {Bkl ■ I < k < I < n} L) {— J n }. Let C(V S ) represent the convex cone associated with the 
basis matrices in V s . Also let C s represent the convex cone of matrices with s being their dominant 
eigenvector (i.e the eigenvector corresponding to the maximal eigenvalue). Then for any R G /C(s), there 
exists ao >0 such that for all a > olq, 

R + ass 11 eC(V s )uC s . (36) 

The proof of Theorem 3 will be presented in several steps (Theorems 4-7 and thereafter). Note that 
we show that (|36l ) can be satisfied even if s is a hyper point of UQP (satisfying d34l). However, since 



March 4, 2013 



DRAFT 



12 



s is the global optimum of UQP for all matrices in C s and C(V a ), the case of ao = can occur only 
when s is a global optimum of UQP associated with R. 

Suppose s is a hyper point of UQP associated with a given positive definite matrix R, and let 6k,i = 
[arg(i?)]fe j j. We define the matrix R + as 

/ |fl(M)|cos(0 fcjJ -(&-&)) 

*MM)=< (37) 
otherwise 

where represents the set of all (fc, Z) such that \9k,i — (4>k — 4>l)\ < 11 1^- Now, let p be a positive real 
number such that 

p> max {\R(k,l)cos(6 kt i - - </>;))!} (38) 

and consider the sequence of matrices {.RW} defined (in an iterative manner) by R^ = R, and 

R (t+i) = R (t) _ {R (t) _ pln><n) Q {ssH) (39) 

for t > 0. The next two theorems (whose proofs are given in the Appendix) study some useful properties 
of the sequence {R^}. 

Theorem 4. {R^} is convergent in at most two iterations: 

R (t) = R (2)^ v t > 2 . (40) 

Theorem 5. Ry^ is a function of p. Let p and p' both satisfy the criterion A38\) . At the convergence of 
{R^} (which is attained for t = 2) we have: 

RW(p') = RW(p) + (p'-p)(ss H ). (41) 

Using the above results, Theorems 6 (whose proof is given in the Appendix) and 7 pave the way for 
a constructive proof of Theorem 3. 

Theorem 6. If s is a hyper point of the UQP associated with R^ = R then it is also a hyper point 
of the UQPs associated with R^ and R^ 2 \ Furthermore, s is an eigenvector of R^ corresponding to 
the eigenvalue np. 

Theorem 7. If s is a hyper point of UQP for R^ = R then it will be the dominant eigenvector of 
R^ if p is sufficiently large. In particular, let p be the largest eigenvalue of R^ which belongs to an 
eigenvector other than s. Then for any p > p/n, s is a dominant eigenvector of R^ 2 \ 
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Proof: We know from Theorem 6 that s is an eigenvector of R^ corresponding to the eigenvalue 
np. However, if s is not the dominant eigenvector of R^ 2 \ Theorem 5 implies that increasing p would 
not change any of the eigenvalues/vectors of R^ except that it increases the eigenvalue corresponding 
to s. As a result, for s to be the dominant eigenvector of R^ we only need p to satisfy np > p or 
equivalently p> p/n, which concludes the proof. ■ 

Returning to Theorem 3, note that R can be written as 

R = R {0) (42) 
= rW + (rW> + R^)Q(ss H )-2pss H . 
For sufficiently large p (satisfying both (I38T ) and the condition of Theorem 7) we have that 

R + 2pss H = R^ + (Rf + flW) ( ss # ) (43) 

where R {2) G C s and (R^ + R^ ] ) Q (ss H ) G C(V S ). Theorem 3 can thus be directly satisfied using 
Eq. (@3]) with a = 2p. 

We conclude this section with two remarks. First of all, the above proof of Theorem 3 does not attempt 
to derive the minimal ao- In the following section we study a computational method to obtain an ao 
which is as small as possible. Secondly, we can use C(V S ) U C s as an approximate characterization of 
JC(s) noting that the accuracy of such a characterization can be measured by the minimal value of ao- 
An explicit formulation of a sub-optimality guarantee for a solution of UQP based on the above /C(s) 
approximation is derived in the following section. 

V. MERIT FOR UQP 

Using the previous results, namely the one-to-one mapping introduced in Theorem 1 and the approx- 
imation of JC(s) derived in Section JV] we build a sequence of matrices (for which the UQP global 
optima are known) whose distance from a given matrix is decreasing. The proposed iterative approach 
can be used to solve for the global optimum of UQP or at least to obtain a local optimum (with an upper 
bound on the sub-optimality of the solution). The sub-optimality guarantees are derived noting that the 
proposed method decreases an upper bound on the sub-optimality of the obtained UQP solution in each 
iteration. 

We know from Theorem 3 that if s is a hyper point of the UQP associated with R then there exist 
matrices Q s G C s , P s G C(V S ) and a scalar «o > such that 

R + a ss H = Q s + P S . (44) 
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Eq. (l44l can be rewritten as 

R + a ss H = (Q 1 + P 1 )Q{ss H ) (45) 

where Q 1 G C\, P\ G C(Vi). We first consider the case of «o = which corresponds to the global 
optimality of s. 

A. Global Optimization of UQP (the Case of oq = 0) 
Consider the optimization problem: 

min \\R-(Q 1 +P 1 )Q(8S H )\\ F (46) 

aeQ n ,Q 1 ee 1 ,Piee(Vi) 

Note that, as C\ U C(V\) is a convex cone, the global optimizers Q x and Pi of (l46l ) for any given s 
can be easily found. On the other hand, the problem of finding an optimal s for fixed R\ = Q 1 + Pi 
is non-convex and hence more difficult to solve globally (see below for details). 

We will assume that R\ is a positive definite matrix. To justify this assumption let R = R (ss H )* 
and note that the eigenvalues of R are exactly the same as those of R, hence R is positive definite. 
Suppose that we have 

/ 

x H Rx > e, V unit-norm x G C nxl 

_ (47) 
\\R-Ri\\f < e 

for some e > 0. It follows from (l47"T ) that 

x Ii R 1 x > x H Rx - \x H Rx - x H R!x\ (48) 

> £ - \x H (R- R±)x\ 

> e-\a\(R-Ri)\ 

> e- \\R-Ri\\ F > 

which implies that R\ is also a positive definite matrix. The conditions in (l47l can be met as follows. 
By considering only the component of R\ in C(V\) (namely Pi) we observe that any positive (i.e. 
with A > 0) diagonal loading of R, which leads to the same diagonal loading of R (as R + XI = 
RQ (ss 11 )* + XI = (R+ XI) (ss H )*), will be absorbed in Pi. Therefore, a positive diagonal loading 
of R does not change \\R — Ri\\f but increases x H Rx by A. We also note that due to \\R — -Ri||_p 
being monotonically decreasing through the iterations of the method, if the conditions in (l47"T ) hold for 
the solution obtained in any iteration, it will hold for all the iterations afterward. 
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In the following, we study a suitable diagonal loading of R that ensures meeting the conditions in 
d4Tb . Next the optimization of the function in (|46"1 ) is discussed through a separate optimization over the 
three variables of the problem. 

• Diagonal loading of R: As will be explained later, we can compute Q x and Pi, (hence R\ = Qx+P\) 
for any initialization of s. In order to guarantee the positive definiteness of R\, define 

e = ||S -Jl! || jr. (49) 

Then we suggest to diagonally load R with A > Ao = —o- n {R) + e$: 

R^ R+XI. (50) 

• Optimization with respect to Q 1 : We restate the objective function of (l46l ) as 

\\R-(Q 1 + P 1 )Q(ss H )\\ F (51) 
= ||( J R0( SS H )*-Pi)-Q 1 || F . 

V v ' 

Rq 

Given Rq, (|46l ) can be written as 

min \\R Q - Qx\\f- (52) 

In EUl . the authors have derived an explicit solution for the optimization problem 

mm\\RQ - QxWf (53) 
Qi 

s.t. Q x l = pi. (p =given) 
The explicit solution of (l53l is given by 

Q x (p) = pI n + (I n -^)(R Q - P I n )(I n -^) (54) 

n n 

p 2 1 

= Rq "I Inxn (-Rolnxn) H 9 (Inxn-Ro Inxn) 

n n n z 

Note that 

Ql(p') " Ql(p) = (P' " P)(lnxl/V^)(lnxl/^) T (55) 

which implies that except for the eigenpair (l n xi/y/n, p), all other eigenvalue/vectors are independent 
of p. Let pq represent the maximal eigenvalue of Qi(0) corresponding to an eigenvector other than 
lnxi/\/n. Therefore, (l52l is equivalent to 

mm \\R Q - (56) 

p 

S.t. p> p . 
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It follows from ([54]) that 

2 



\Rq-Qi(p)\\ 2 f = J2' 



n 

k=i 

Jh 



p 2G k H 

n n n 2 



(57) 



where Gk and H are the sum of the k row and, respectively, the sum of all entries of Rq. The p that 
minimizes (1571 ) is given by 



P 

n 



l^(^k--)=- (58) 

fe=i v 7 



which implies that the minimizer p = p* of (1561) is equal to 



— — > Po, 



P*={ n n - ' (59) 
po otherwise. 



Finally, the optimal solution Q x to (1521 ) is given by 

Qi = Qxip*)- (60) 

• Optimization with respect to P\: Similar to the previous case, (|46l can be rephrased as 



min \\Rp-PAf (61) 
Q 1 eC(V 1 ) 

where Rp = RQ (ss H )* — Q 1 . The solution of doTt is simply given by 

| R' P (k,l) R' P (k,l)>0ork = l, 
Pi{k,l) = < (62) 
otherwise 

where R' P = ^{R P }. 

• Optimization with respect to s: Suppose that Q x and P\ are given and that R\ = Q t + P\ is a 
positive definite matrix (see the discussion on this aspect following Eq. (l46ll). We consider a relaxed 
version of (|46l , 

min (aiaf )\\p (63) 

The objective function in (l63l can be re-written as 

|| J R-.R 1 Q(s 1 sf)||! (64) 
= Hfl-Diag^OiiiDiag^lll 

= tr^ 2 ) + tr(^) - 2&{tr(.R Diag(si) #i Diag(s^))}. 
Note that only the third term of (|64l is a function of si and S2- Moreover, it can be verified that ET1 

tr(R Diag(si) i?i Diag(s^)) = s$(R © R% )ai- (65) 
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TABLE I 
The MERIT Algorithm 



(A) The case of Qo = 

Step 0: Initialize the variables Q 1 and Pi with J. Let s be a random vector 

in n n . 

Step 1: Perform the diagonal loading of R as in d49t-(150t (note that this 
diagonal loading is sufficient to keep Pi = Q x + Pi always positive definite). 
Step 2: Obtain the minimum of J46t with respect to Q 1 as in d60t . 
Step 3: Obtain the minimum of l !46t with respect to Pi using l !62b . 
Step 4: Minimize l |46l > with respect to s using l |66t . 

Step 5: Goto step 2 until a stop criterion is satisfied, e.g. \\R — (Q 1 + Pi) 
{ss H ) || f < eo (or if the number of iterations exceeded a predefined maximum 
number). 

(B) The case of cto > 

Step 0: Initialize the variables (s, Q ± , P±) using the results obtained by the 
optimization of d46b as in Table HJ-A. 

Step 1: Set 8 (the step size for increasing ao in each iteration). Let So be the 
minimal 8 to be considered and ao = 0. 

Step 2: Let ag" = a , aT w = Qo + 8 and R' = R + a r lew ss H . 

Step 3: Solve d67t using the steps 2-5 in Table [JA (particularly step 4 must 

be applied to d69|i). 

Step 4: If \\R' - (Q ± + Pi) (ss H )\\ F < e do: 

• Step 4-1: If 8 > So, let 8 8/2 and initialize ( 167b with the previously 
obtained variables (s, Q 1; Pi) for Qo = ci^ re . Goto step 2. 

• Step 4-2: If S < 5 , stop. 
Else, let ao = a lew and goto step 2. 



As R is positive definite, we can employ the power method-like iterations introduced in (l32l i to 
obtain a solution to d46l ) i.e. starting from the current s = s' ', a local optimum of the problem can be 
obtained by the iterations 

s (t+l) = e j a-rgdReR^sW) _ (66) 

Finally, the proposed algorithmic optimization of (l46l i based on the above results is summarized in 
Table UA. 
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B. Achieving a Local Optimum of UQP (the Case of a® > 0) 

There exist examples for which the objective function in (l46"T ) does not converge to zero. As a result, 
the proposed method cannot obtain a global optimum of UQP in such cases. However, it is still possible 
to obtain a local optimum of UQP for some ao > 0. To do so, we solve the optimization problem, 

min \\R' - (Q x + Pi) (ss H )\\ F (67) 

sen,Q 1 eCi,.Pi€C(Vi) 

with R = R + aoss 11 , for increasing ao- The above optimization problem can be tackled using the 
same tools as proposed for (l46l . In particular, note that increasing ao decreases (I67T ). To observe this, 
suppose that the solution (s,Qi,Pi) of (l67l ) is given for an ao > 0. The minimization of (167T ) with 
respect to Q 1 for a.Q ew = a + 5 {5 > 0) yields Q 1 G Ci such that 

\\R + a% ew ss H - (Q 1 + Pi) (ss H )\\ F (68) 

< \\R + a^ w 88 H - ((Q 1 + 511 T ) + Pi) (ss H )\\ F 

= \\R + a ss H -(Q 1 + P 1 )Q{ss H )\\ F 

where Q x + 511 T G Ci. The optimization of (I67T ) with respect to Pi can be dealt with as before (see 
(|46j ) and it leads to a further decrease of the objective function. Furthermore, 

\\R + a Q 88 H -(Q 1 + P 1 )Q(88 H )\\ F (69) 

= ||P + X'l - (Qi + Pi - a ll T + X'l) (ss H )\\ F 

which implies that a solution s of (l67l ) can be obtained via optimizing (|69l with respect to s in a similar 
way as we described for (|46l provided that A' > is such that Q x + Pi — aoll T + X'l is positive definite. 
Finally, note that the obtained solution (s,Qi,Pi) of (l46l ) can be used to initialize the corresponding 
variables in (I67T ). In effect, the solution of (l67l ) for any ao can be used for the initialization of (l67l i with 
an increased ao- 

Based on the above discussion and the fact that small values of ao are of interest, a bisection approach 
can be used to obtain ao- The proposed method for obtaining a local optimum of UQP along with the 
corresponding ao is described in Table [Q-B. 

C. Sub-optimality Analysis 

In this sub-section, we show that the proposed method can provide a sub-optimality guarantee (7) that 
is close to 1. Let ao = (as a result R = R) and define 

P = P' - (Qi + Pi) (ss 11 ) (70) 
> „ ' 

R. 
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where Q t G C\ and Pi G C(Vi). The global optimum of the UQP associated with R s is s. We have 
that 

max s' H Rs' < max s' H R s s' + max s' H Es' (71) 

s'esi" s'esi" s'esi" 

< max s' H R s s' + naAE) 
= s H R s s + nai(E) 



Furthermore, 



max s' H Rs' > max s' n R s s' + min rEs' (72) 

s'ef!" s'efi n s'efi n 



> max s' H R s s' + na n (E) 



= s^.R s s + racr„( J E) 

As a result, an upper bound and a lower bound on the objective function for the global optimum of (l46l ) 
can be obtained at each iteration. Furthermore, as 

Wx{E)\ < \\E\\ F , \a n (E)\ < \\E\\ F (73) 

if \\E\\ F converges to zero we conclude for (1711 and (l72l that 



max s' H Rs' = s H R s s = s H Rs (74) 

and hence s is the global optimum of the UQP associated with R (i.e. a sub-optimality guarantee of 
7 = 1 is achieved). 

Next, suppose that we have to increase «o i n order to obtain the convergence of to zero. In such a 

case, we have that R = R s — aoss H and as a result, max s / e f^ s' H R s s' — a^n 2 < max s / e f^ s' H Rs' < 
max s / e si» s' H R s s' or equivalently, 

s H R s s - a n 2 < max s' H Rs' < s H R s s. (75) 



The provided sub-optimality guarantee is thus given by 

s Rs aon 2 
s H R*s s H R*s 



^ = -JT5- = 1 -7H7rz- (76) 



Note that while solving the optimization problem (1671 ) does not necessarily yield the exact optimal 
solution to UQP, the so-obtained solution can be still optimal. We also note that (l76l ) generally yields 
tighter sub-optimality guarantees than the currently known approximation guarantee (i.e. 7r/4 for SDR). 
The following section provides empirical evidence for such a fact. 
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Fig. 1. Different metrics versus the iteration number for an UQP solved by MERIT, (a) the UQP objective corresponding to the 
true matrix (it), the approximated matrix (it s ) and also the upper/lower bounds at each iteration. The sub-optimality bounds 
are updated using ilH-il2). (b) the criterion = \\R — R 3 \\f (it reaches values which are practically zero). 



VI. Numerical Examples 

In order to examine the performance of the proposed method, several numerical examples will be 
presented. Random Hermitian matrices R are generated using the formula 

n 

# = J> fc af (77) 
fc=i 

where {x^} are random vectors in C n whose real-part and imaginary-part elements are i.i.d. with a 
standard Gaussian distribution A/"(0, 1). In all cases, we stopped the iterations when ||.E||f < 10 -9 . 

We use the MERIT algorithm to solve the UQP for a random positive definite matrix of size n = 16. 
The obtained values of the UQP objective for the true matrix (R) and the approximated matrix (R s ) as 
well as the sub-optimality bounds (derived in (TTTb and (|72T i) are depicted in Fig. 1 versus the iteration 
number. In this example, a sub-optimality guarantee of 7 = 1 is achieved which implies that the method 
has successfully obtained the global optimum of the considered UQP. A computational time of 3.653 sec 
was required on a standard PC to accomplish the task. 

Next, we solve the UQP for 20 full-rank random positive definite matrices of sizes n £ {8, 16, 32, 64}. 
Inspired by ifTTI and 11221 . we also consider rank-deficient matrices R = J2k=i x k x k where {xk} are as 
in (1771 ). but d <C n. The performance of MERIT for different values of d is shown in Table [TIJ Interestingly, 
the solution of UQP for rank-deficient matrices appears to be obtained more efficiently than for the full- 
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n 


Rank (d) 


#problems for 


Average 7 


Minimum 7 


Average CPU 


#problems 






which 7=1 






time (sec) 


solved by SDR 


8 


2 


17 


0.9841 


0.8184 


0.13 


4 




8 


16 


0.9912 


0.9117 


0.69 


7 




2 


15 


0.9789 


0.8301 


1.06 


2 


16 


4 


13 


0.9773 


0.8692 


1.58 


10 




16 


4 


0.9610 


0.8693 


3.54 


13 




2 


g 


9536 


0.8190 


47.04 


3 


32 


6 


4 


0.9077 


0.8106 


55.59 


7 




32 


2 


0.9031 


0.8021 


94.90 


16 




2 


3 


0.8893 


0.8177 


406.56 


4 


64 


8 


1 


0.8567 


0.7727 


560.35 


10 




64 





0.8369 


0.7811 


1017.69 


15 



TABLE II 

Comparison of the performance of MERIT (see Table[U) and SDR ifTTll when solving the UQP for 20 random 

POSITIVE DEFINITE MATRICES OF DIFFERENT SIZES n AND RANKS d. 



rank matrices. For each problem solved by MERIT, we also let the SDR algorithm of ifTTl use the same 
computational time for solving the problem. The SDR algorithm is able to solve the problem only if its 
core semi-definite program can be solved within the available time. Any remaining time is dedicated to 
the randomization procedure. The results can be found in TablelTTl Note that the maximum UQP objective 
values obtained by MERIT and SDR were nearly identical in those cases in which SDR was able to 
solve the UQPs in the same amount of time as MERIT. Note also that given the solutions obtained by 
MERIT and SDR as well as the sub-optimality guarantee of MERIT, a case-dependent sub-optimality 
guarantee for SDR can be computed as 

A f VSDR \ 

7SDR = 7MERIT ) ■ ™ 

This can be used to examine the goodness of the solutions obtained by SDR. 

Besides random matrices, we also consider several other matrix structures for which solving the UQP 
using the proposed method is not "hard", as explained below. 

• Case 1: An exponentially shaped disturbance matrix fll] with correlation coefficient rj = 0.8, 

M(k,l) =V fc-{| , l<k,l<n. (79) 
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Case 2: A disturbance matrix with the structure 

M(k, I) = r)[ k - ll e j2wp{k - l) + Wi] l 2 k ~ 11 + 1(T 2 I(£:, I), 1 < k, I < n (80) 

whose terms represent the effects of sea clutter, land clutter and thermal noise, respectively. The 

values of (771, 772, p) are set to (0.8,0.9,0.2) in accordance to an example provided in ll23l . 

Case 3: A disturbance matrix accounting for both discrete clutter scatterers and thermal noise lUTI , 



M = Y,^Pv^Pv dtk +r}I (81) 

k=l 

where n c = 10, % = 10 3 , v d)k = (k - l)/2, 

p Vdk = il.r'"-' , e iM"-lK, fc) T ; 1 < k < ncj (g2) 

and r\ = 10~ 2 . The chosen values are the same as those considered in ifTTl . 
We let R = M^ 1 (pp H )* (see ([3]) and the following discussion) where p is an unimodular vector 
with a structure similar to that of {p Vd k } in (f82l . The UQP for the above cases is solved via MERIT 
using 20 different random initializations for sizes n € {8, 16, 32, 64}. Similar to the previous example, 
we also used SDR to solve the same UQPs. The results are shown in Table Jill The obtained solutions 
can be considered to be quite accurate in the sense of a sub-optimality guarantee 7 close to one. 

A different code design problem arises when synthesizing waveforms that have good resolution prop- 
erties in range and Doppler [3]-[5],[24]-||26]. In the following, we consider the design of a thumbtack 
CAF (see the definitions in sub-section II-AI ): 

otherwise. 

Suppose n = 53, let T be the time duration of the total waveform, and let t p = T/n represent the time 
duration of each sub-pulse. Define the weighting function as 

w(t, f) = < (84) 
otherwise, 

where * = [-10t p , 10i p ]x [-2/T, 2/T] is the region of interest and ^ mt = ({-t p ,t p ]\{0})x({-l/T,l/T} 
\{0}) is the mainlobe area which is excluded due to the sharp changes near the origin of d(r, /). Note 
that the time delay r and the Doppler frequency / are typically normalized by T and 1/T, respectively, 
and as a result the value of t p can be chosen freely without changing the performance of CAF design. 
The synthesis of the desired CAF is accomplished via the cyclic minimization of © with respect to x 
and y (see sub-section II-AI ). In particular, we use MERIT to obtain a unimodular x in each iteration. A 
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TL 


Rank (d) 


#problems for 


Average *y 


Minimum *y 


Average CPU 


#problems 






which 7=1 






time (sec) 


solved by SDR 




Cg.sc 1 


20 


1.0000 


1.0000 


2.82 


17 


8 


Case 2 


20 


1.0000 


1.0000 


0.60 


20 




Case 3 


20 


1.0000 


1.0000 


0.27 


10 




\^ase j. 


20 


1.0000 


1.0000 


42.83 


20 


16 


Case 2 


18 


0.9812 


0.8075 


21.58 


20 




Case 3 


20 


1.0000 


1.0000 


2.01 


12 




r 

y^ase j 


20 


1.0000 


1.0000 


990 90 


20 


32 


Case 2 


19 


0.9995 


0.9913 


525.34 


20 




Case 3 


20 


1.0000 


1.0000 


7.52 


7 




Case 1 


17 


0.9901 


0.9862 


5574.98 


20 


64 


Case 2 


16 


0.9540 


0.8359 


2053.26 


20 




Case 3 


20 


1.0000 


1.0000 


22.78 


9 



TABLE III 

Comparison of the performance of MERIT (see TableQ} and SDR fTTll when solving the UQP for the 

MATRIX STRUCTURES DESCRIBED IN CASES 1-3 USING 20 DIFFERENT INITIALIZATIONS AND FOR DIFFERENT SIZES n. 



Bjorck code is used to initialize both vectors x and y. The Bjorck code of length n = p (where p is a 
prime number for which p = l (mod 4)) is given by b(k) = e j{ p ] arccos ( 1 /(H- v / p)) ) o < k < p, with (|) 
denoting the Legendre symbol. Fig. 2 depicts the normalized CAF modulus of the Bjorck code (i.e. the 
initial CAF) and the obtained CAF using the UQP formulation in (TTTb and the proposed method. Despite 
the fact that designing CAF with a unimodular transmit vector x is a rather constrained problem, MERIT 
is able to efficiently suppress the CAF sidelobes in the region of interest. 

VII. Concluding Remarks 

A computational approach to the NP-hard problem of optimizing a quadratic form over the unimodular 
vector set (called UQP) has been introduced. The main results can be summarized as follows: 

• Some applications of the UQP were reviewed. It was shown that the solution to UQP is not 
necessarily unique. Several examples were provided for which an accurate global optimum of UQP 
can be obtained efficiently. 

• Using a relaxed version of UQP, a specialized local optimization scheme for UQP was devised and 
was shown to yield superior results compared to any general local optimization of UQP. 
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\x(tJ)\ | X (r,/)j (dB) 




(b) left: 3D plot, right: 2D plot 

Fig. 2. The normalized CAF modulus for (a) the Bjorck code of length n — 53 (i.e. the initial CAF), and (b) the UQP 
formulation in (TT) and MERIT. 

• It was shown that the set of matrices (X(s)) leading to the same solution (s) as the global optimum of 
UQP is a convex cone. An one-to-one mapping between any two such convex cones was introduced 
and an approximate characterization of JC(s) was proposed. 

• Using the approximate characterization of IC(s), an iterative approach (called MERIT) to the UQP 
was proposed. It was shown that MERIT provides case-dependent sub-optimality guarantees. The 
available numerical evidence shows that the sub-optimality guarantees obtained by MERIT are 
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generally better than the currently known approximation guarantee (of 7r/4 for SDR). 
• Numerical examples were provided to examine the potential of MERIT for different UQPs. In 
particular, it was shown that the UQP solutions for certain matrices used in active sensing code 
design can be obtained efficiently via MERIT. 
We should note that no theoretical efficiency assessment of the method was provided. It is clear that 
C(V a ) U C s C JC(s). A possible approach would be to determine how large is the part of IC(s) that is 
"covered" by C(V a ) UC S . However, this problem is left for future work. Furthermore, a study of m-UQP 
using the ideas in this paper will be the subject of another paper. 



A. Proof of Theorem 1 

In order to verify the first part of the theorem, consider any two matrices R\,R2 G JC(s). For any 
two non-negative scalars 71 , 72 we have that 



Clearly, if some s = s is the global maximizer of both s R\s and s J?2« then it is the global maximizer 
of s H (jiR\ + j 2 R 2 )s which implies 71-R1 + 72-R2 G JC(s). 
The second part of the theorem can be shown noting that 



= sfRsi 

for all si,s 2 G Vt n and s = s* s 2 . Therefore, if R G K,(s\) then R © (sqSq ) G /C(s 2 ) (for 
So = si 82) and vice versa. 

B. Proof of Theorem 2 

It is well-known that x H Rx < cillaslln for all vectors x G C n . Let 



Appendix 



s H (j 1 R 1 + 72-R 2 )s = - fl s H R 1 s + -f 2 s H R 2 s. 



(85) 



sf (fl©(s *o ))*2 



(Sq S2) R(Sq 82) 



(86) 




(87) 



It follows from (1241 ) that s = Ua' and therefore 



s H Rs = a' H Y±a' = aAa! 



(88) 



0l||a||! = no ~i 



which implies the global optimality of s for the considered UQP. 
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C. Proof of Theorem 4 

It is worthwhile to observe that the convergence rate of {R^} is not dependent on the problem 
dimension (n), as each entry of {R^} is treated independently from the other entries (i.e. all the 
operations are element-wise). Therefore, without loss of generality we study the convergence of one 
entry (say {R^\k, I)} = {r t e jdt }) in the following. 

Note that in cases for which \9 t — (4>k — 4>i)\ > 7r /2, the next element of the sequence \rte? 6t } can be 
written as 

r t+1 e jdt+1 = r t e jdt + pe j ^ k ~^ (89) 

which implies that the proposed operation tends to make 9t closer to (fa — 4>i) i n eacn iteration, and 
finally puts 6 t within the tt/2 distance from (fa — fa). 

Let us suppose that \9q — (fa — 4>i)\ > tt/2, and that the latter phase criterion remains satisfied for all 
9 t , t <T. We have that 

r T e jdT = r e je ° + TpeP^"'^ (90) 

which yields 

r T cos(6» T - (fa - fa)) = r cos(9 - (fa - 4>i)) + Tp. (91) 

Therefore it takes only T = \— rocos(#o — (<Pk — fa))/p\ = 1 iteration for 9 t to stand within the tt/2 
distance from (fa — fa). 

Now, suppose that \9q — (fa^ — 4>i)\ < tt/2. For every t > 1 we can write that 

r t+1 e jf>t+1 = r t e jdt + pe j ^ k ~ M (92) 

- r t cos(e t -(fa-fa))ei^-^ 

Let 5t+i = r t +ie^ t+1 — r t e^ et . The first equality in (|92l implies that 

St+i = e^ k -^(p - r t cos(9 t - (<f> k - 0,)))- (93) 
On the other hand, the second equality in (|92l implies that 

<5t+i = je^-^in sm(9 t - (fa - fa)) (94) 
-r t -i sin(0 t _i - (fa - 4>i))) 
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for all t > 1. Note that in d93l ) and d94l ), is a complex number having different phases. We conclude 

St+i = 0, Vt>l (95) 

which shows that the sequence {rte^ et } is convergent in one iteration. In sum, every entry of the matrix 
R will converge in at most two iterations (i.e. at most one to achieve a phase value within the ir/2 
distance from (4>k — 4>i), and one iteration thereafter). 

D. Proof of Theorem 5 

We use the same notations as in the proof of Theorem 4. If \9q — (<pk — <j>i)\ < tt/2 then 

r 2 e j92 = ne jei (96) 
= r e? e ° + peP^"-^ 

- roco S (0o-(0 fc -^))e i(0fc -* ) . 

On the other hand, if \Q - {<f> k - <fii)\ > ir/2 we have that rie jdl = r e je " + j oe^ h-< ^. As a result, 
ri cos(0i - (<f> k - fa)) = p + r cos(6> - (0& - <j>i)) which implies 

r 2 e jd2 = r ie jdl + p e ^ k ~^ (97) 

- ncos(ei-(^-^))e»'(**-^ 
= r e je " + pe j ^ k ^ l) 

- r QC os{e Q -{4> k -^))e^^\ 
Now, it is easy to verify that (|4TT > follows directly from (|96l and 



Proof of Theorem 6 

If s is a hyper point of UQP associated with R^ = R then we have that arg(s) = axg(Rs). Let 
Rs = v s where v is a non-negative real-valued vector in R n . It follows that 

n 

v(k)e j ^ k = \R{k,l)\e jdk - l e j<t>l (98) 
i=i 

or equivalently 

n 

v(k) = \ R (k, l)^ 6 "- 1 -^"-^ (99) 
i=i 
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which implies that 

' Er=ii«(^oicos(^-(^-^))>o 

_ EF=i 01 M0k,i - {<t>k - </>/)) = o 

for all I < k < n. Now, note that the recursive formula of the sequence {R^'} can be rewritten as 

R (t+i) = R (t) _ Diag(s ) (jR W _ plnxn) Diag( S *) (101) 

and as a result, 

R (t+i) s = R (t) s _ Diag(s ) (Ji W _ p i nxn ) i nxl . (102) 

It follows from (11021 ) that if s is a hyper point of the UQP associated with R^ (which implies the 
existence of non-negative real- valued vector such that R^s = v^Qs), then there exists G W 1 

for which = © s and therefore, 

n 

v^ik)^^ = ^2\R^(k,l)\eP e "' l eP' t>l (103) 



l=i 

/ / n 
J=l 



np e J0 \ 



Eq. (11031 ) can be rewritten as 

n 

V ^ t+1 \k) = ^|jtf t )(M)|e' ( * w ~ ( *" , ~* ,)) (104) 



i=i 

/ n 



As indicated earlier, s being a hyper point for R^ assures that the imaginary part of (11041 ) is zero. To 
show that s is a hyper point of the UQP associated with R^ t+1 \ we only need to verify that i>( t+1 ) [k) > 0: 

n 

v {t+1) (k) = J2\R^(k,l)\cae(9 k ,i-(<f> k -<l>i)) (105) 
i=i 

= np 

+ Yl l^ (i) (^,0|cos(^-(^-^)) 

I: (k,l)<£e 

Now note that the positivity of i>( <+1 )(/c) is concluded from (|38l l. In particular, based on the discussions 
in the proof of Theorem 4, for t = 1, there is no such that \6 k i — {4> k — (pi)\ > 7r/2 and therefore 
i/ 2 ) (&;) = for all 1 < k < n. As a result, 

R^s = nps (106) 
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which implies that s is an eigenvector of corresponding to the eigenvalue np. 
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